1: Ordination

#Load Packages and Prepare Data

library("phyloseq"); packageVersion("phyloseq")
## [1] '1.32.0'
data(GlobalPatterns)
library("ggplot2"); packageVersion("ggplot2")
## [1] '3.3.2'
library("plyr")

theme_set(theme_bw())

Filter low-occurence OTUs

GP = GlobalPatterns
wh0 = genefilter_sample(GP, filterfun_sample(function(x) x > 5), A=0.5*nsamples(GP))
GP1 = prune_taxa(wh0, GP)

Even sampling Depth

GP1 = transform_sample_counts(GP1, function(x) 1E6 * x/sum(x))

Only Keep Five Mos Abundant Phyla

phylum.sum = tapply(taxa_sums(GP1), tax_table(GP1)[, "Phylum"], sum, na.rm=TRUE)
top5phyla = names(sort(phylum.sum, TRUE))[1:5]
GP1 = prune_taxa((tax_table(GP1)[, "Phylum"] %in% top5phyla), GP1)

Create human vs nonhuman categorical variable

human = get_variable(GP1, "SampleType") %in% c("Feces", "Mock", "Skin", "Tongue")
sample_data(GP1)$human <- factor(human)

Four Main ORdination Plots

1) Just OTUs

GP.ord <- ordinate(GP1, "NMDS", "bray")
## Square root transformation
## Wisconsin double standardization
## Run 0 stress 0.1333468 
## Run 1 stress 0.1648238 
## Run 2 stress 0.1577391 
## Run 3 stress 0.1333468 
## ... Procrustes: rmse 3.039838e-06  max resid 7.386318e-06 
## ... Similar to previous best
## Run 4 stress 0.161634 
## Run 5 stress 0.2544504 
## Run 6 stress 0.1333468 
## ... Procrustes: rmse 1.045347e-06  max resid 1.726379e-06 
## ... Similar to previous best
## Run 7 stress 0.157644 
## Run 8 stress 0.1748452 
## Run 9 stress 0.153803 
## Run 10 stress 0.2838378 
## Run 11 stress 0.1333468 
## ... Procrustes: rmse 1.722968e-05  max resid 3.041404e-05 
## ... Similar to previous best
## Run 12 stress 0.1622613 
## Run 13 stress 0.1333468 
## ... Procrustes: rmse 2.197773e-06  max resid 6.541958e-06 
## ... Similar to previous best
## Run 14 stress 0.1333468 
## ... New best solution
## ... Procrustes: rmse 4.018371e-06  max resid 1.157531e-05 
## ... Similar to previous best
## Run 15 stress 0.1674589 
## Run 16 stress 0.1333468 
## ... Procrustes: rmse 5.693894e-06  max resid 1.688925e-05 
## ... Similar to previous best
## Run 17 stress 0.1333468 
## ... Procrustes: rmse 3.422433e-06  max resid 9.37327e-06 
## ... Similar to previous best
## Run 18 stress 0.146027 
## Run 19 stress 0.1529609 
## Run 20 stress 0.1625845 
## *** Solution reached
p1 = plot_ordination(GP1, GP.ord, type="taxa", color="Phylum", title="taxa")
print(p1)

Simplify Visually

p1 + facet_wrap(~Phylum, 3)

2) Just Samples

p2 = plot_ordination(GP1, GP.ord, type="samples", color="SampleType", shape="human") 
p2 + geom_polygon(aes(fill=SampleType)) + geom_point(size=5) + ggtitle("samples")

3) Biplot Graphic

p3 = plot_ordination(GP1, GP.ord, type="biplot", color="SampleType", shape="Phylum", title="biplot")
# Some stuff to modify the automatic shape scale
GP1.shape.names = get_taxa_unique(GP1, "Phylum")
GP1.shape <- 15:(15 + length(GP1.shape.names) - 1)
names(GP1.shape) <- GP1.shape.names
GP1.shape["samples"] <- 16
p3 + scale_shape_manual(values=GP1.shape)

4) Split Graphic

p4 = plot_ordination(GP1, GP.ord, type="split", color="Phylum", shape="human", label="SampleType", title="split") 
p4

Set to Default Color Scale

gg_color_hue <- function(n){
    hues = seq(15, 375, length=n+1)
    hcl(h=hues, l=65, c=100)[1:n]
}
color.names <- levels(p4$data$Phylum)
p4cols <- gg_color_hue(length(color.names))
names(p4cols) <- color.names
p4cols["samples"] <- "black"
p4 + scale_color_manual(values=p4cols)

Supported Ordination Methods

dist = "bray"
ord_meths = c("DCA", "CCA", "RDA", "DPCoA", "NMDS", "MDS", "PCoA")
plist = lapply(as.list(ord_meths), function(i, physeq, dist){
        ordi = ordinate(physeq, method=i, distance=dist)
        plot_ordination(physeq, ordi, "samples", color="SampleType")
}, GP1, dist)
## Square root transformation
## Wisconsin double standardization
## Run 0 stress 0.1333468 
## Run 1 stress 0.1385323 
## Run 2 stress 0.156695 
## Run 3 stress 0.1333468 
## ... New best solution
## ... Procrustes: rmse 1.296711e-06  max resid 3.456084e-06 
## ... Similar to previous best
## Run 4 stress 0.1617642 
## Run 5 stress 0.1626327 
## Run 6 stress 0.1333468 
## ... Procrustes: rmse 1.762521e-06  max resid 4.140732e-06 
## ... Similar to previous best
## Run 7 stress 0.1522771 
## Run 8 stress 0.1529686 
## Run 9 stress 0.146027 
## Run 10 stress 0.1820691 
## Run 11 stress 0.1569043 
## Run 12 stress 0.1625845 
## Run 13 stress 0.146027 
## Run 14 stress 0.1333468 
## ... Procrustes: rmse 3.472052e-06  max resid 1.03564e-05 
## ... Similar to previous best
## Run 15 stress 0.1529686 
## Run 16 stress 0.3021346 
## Run 17 stress 0.1449078 
## Run 18 stress 0.17179 
## Run 19 stress 0.1333468 
## ... Procrustes: rmse 8.106219e-06  max resid 2.431961e-05 
## ... Similar to previous best
## Run 20 stress 0.1333468 
## ... New best solution
## ... Procrustes: rmse 1.75444e-06  max resid 4.440386e-06 
## ... Similar to previous best
## *** Solution reached
names(plist) <- ord_meths
pdataframe = ldply(plist, function(x){
    df = x$data[, 1:2]
    colnames(df) = c("Axis_1", "Axis_2")
    return(cbind(df, x$data))
})
names(pdataframe)[1] = "method"
p = ggplot(pdataframe, aes(Axis_1, Axis_2, color=SampleType, shape=human, fill=SampleType))
p = p + geom_point(size=4) + geom_polygon()
p = p + facet_wrap(~method, scales="free")
p = p + scale_fill_brewer(type="qual", palette="Set1")
p = p + scale_colour_brewer(type="qual", palette="Set1")
p

plist[[2]] 

p = plist[[2]] + scale_colour_brewer(type="qual", palette="Set1")
p = p + scale_fill_brewer(type="qual", palette="Set1")
p = p + geom_point(size=5) + geom_polygon(aes(fill=SampleType))
p

MDS (“PCoA”) on Unifrac Distances

ordu = ordinate(GP1, "PCoA", "unifrac", weighted=TRUE)
plot_ordination(GP1, ordu, color="SampleType", shape="human")

p = plot_ordination(GP1, ordu, color="SampleType", shape="human")
p = p + geom_point(size=7, alpha=0.75)
p = p + scale_colour_brewer(type="qual", palette="Set1")
p + ggtitle("MDS/PCoA on weighted-UniFrac distance, GlobalPatterns")

2: Alpha Diversity Graphics

Load packages, set parameters

library("phyloseq"); packageVersion("phyloseq")
## [1] '1.32.0'
data("GlobalPatterns")
library("ggplot2"); packageVersion("ggplot2")
## [1] '3.3.2'
theme_set(theme_bw())
pal = "Set1"
scale_colour_discrete <-  function(palname=pal, ...){
  scale_colour_brewer(palette=palname, ...)
}
scale_fill_discrete <-  function(palname=pal, ...){
  scale_fill_brewer(palette=palname, ...)
}

Prepare Data

GP <- prune_species(speciesSums(GlobalPatterns) > 0, GlobalPatterns)
## Warning: 'prune_species' is deprecated.
## Use 'prune_taxa' instead.
## See help("Deprecated") and help("phyloseq-deprecated").
## Warning: 'speciesSums' is deprecated.
## Use 'taxa_sums' instead.
## See help("Deprecated") and help("phyloseq-deprecated").

Plot Examples

plot_richness(GP)

plot_richness(GP, measures=c("Chao1", "Shannon"))

plot_richness(GP, x="SampleType", measures=c("Chao1", "Shannon"))

sampleData(GP)$human <- getVariable(GP, "SampleType") %in% c("Feces", "Mock", "Skin", "Tongue")
## Warning: 'getVariable' is deprecated.
## Use 'get_variable' instead.
## See help("Deprecated") and help("phyloseq-deprecated").
## Warning: 'sampleData' is deprecated.
## Use 'sample_data' instead.
## See help("Deprecated") and help("phyloseq-deprecated").
## Warning: 'sampleData<-' is deprecated.
## Use 'sample_data<-' instead.
## See help("Deprecated") and help("phyloseq-deprecated").
plot_richness(GP, x="human", color="SampleType", measures=c("Chao1", "Shannon"))

GPst = merge_samples(GP, "SampleType")
# repair variables that were damaged during merge (coerced to numeric)
sample_data(GPst)$SampleType <- factor(sample_names(GPst))
sample_data(GPst)$human <- as.logical(sample_data(GPst)$human)
p = plot_richness(GPst, x="human", color="SampleType", measures=c("Chao1", "Shannon"))
p + geom_point(size=5, alpha=0.7)

More about GGPLOT2

Check the lists present

p$layers
## [[1]]
## geom_point: na.rm = TRUE
## stat_identity: na.rm = TRUE
## position_identity 
## 
## [[2]]
## mapping: ymax = ~value + se, ymin = ~value - se 
## geom_errorbar: na.rm = FALSE, orientation = NA, width = 0.1, width = 0.1, flipped_aes = FALSE
## stat_identity: na.rm = FALSE
## position_identity
p$layers <- p$layers[-1]
p + geom_point(size=5, alpha=0.7)

3: Heatmap Plots

Examples

library("phyloseq"); packageVersion("phyloseq")
## [1] '1.32.0'
library("ggplot2"); packageVersion("ggplot2")
## [1] '3.3.2'
theme_set(theme_bw())

Plot a 300-Taxa Dataset

data("GlobalPatterns")
gpt <- subset_taxa(GlobalPatterns, Kingdom=="Bacteria")
gpt <- prune_taxa(names(sort(taxa_sums(gpt),TRUE)[1:300]), gpt)
plot_heatmap(gpt, sample.label="SampleType")
## Warning: Transformation introduced infinite values in discrete y-axis

Subset a smaller dataset based on an Archaeal phylum

gpac <- subset_taxa(GlobalPatterns, Phylum=="Crenarchaeota")

Default plot_heatmap settings

plot_heatmap(gpac)
## Warning: Transformation introduced infinite values in discrete y-axis

Re-label by a sample variable and taxonomic family

(p <- plot_heatmap(gpac, "NMDS", "bray", "SampleType", "Family"))
## Warning: Transformation introduced infinite values in discrete y-axis

Relabel Axis Titles

p$scales$scales[[1]]$name <- "My X-Axis"
p$scales$scales[[2]]$name <- "My Y-Axis"
print(p)
## Warning: Transformation introduced infinite values in discrete y-axis

Now repeat the plot, but change the color scheme

plot_heatmap(gpac, "NMDS", "bray", "SampleType", "Family", low="#000033", high="#CCFF66")
## Warning: Transformation introduced infinite values in discrete y-axis

plot_heatmap(gpac, "NMDS", "bray", "SampleType", "Family", low="#000033", high="#FF3300")
## Warning: Transformation introduced infinite values in discrete y-axis

plot_heatmap(gpac, "NMDS", "bray", "SampleType", "Family", low="#000033", high="#66CCFF")
## Warning: Transformation introduced infinite values in discrete y-axis

Dark on White Color Scheme

plot_heatmap(gpac, "NMDS", "bray", "SampleType", "Family", low="#66CCFF", high="#000033", na.value="white")
## Warning: Transformation introduced infinite values in discrete y-axis

Cream Colored

plot_heatmap(gpac, "NMDS", "bray", "SampleType", "Family", low="#FFFFCC", high="#000033", na.value="white")
## Warning: Transformation introduced infinite values in discrete y-axis

Different Ordination Methods, Distances

Jaccard Distance

plot_heatmap(gpac, "NMDS", "jaccard")
## Warning: Transformation introduced infinite values in discrete y-axis

Detrended Correspodence Analysis

plot_heatmap(gpac, "DCA", "none", "SampleType", "Family")
## Warning: Transformation introduced infinite values in discrete y-axis

Unconstrained Reduncancy Analysis

plot_heatmap(gpac, "RDA", "none", "SampleType", "Family")
## Warning: Transformation introduced infinite values in discrete y-axis

PCoA/MDS of Bray-Curtis Distance

plot_heatmap(gpac, "PCoA", "bray", "SampleType", "Family")
## Warning: Transformation introduced infinite values in discrete y-axis

PCoA/MDS of Unifrac Distance

plot_heatmap(gpac, "PCoA", "unifrac", "SampleType", "Family")
## Warning: Transformation introduced infinite values in discrete y-axis

Weighted-UniFrac distance and MDS/PCoA ordination

plot_heatmap(gpac, "MDS", "unifrac", "SampleType", "Family", weighted=TRUE)
## Warning: Transformation introduced infinite values in discrete y-axis

Hierarchichal Clustering

heatmap(otu_table(gpac))

4: Networks

Load Package and Example Data

library(phyloseq); packageVersion("phyloseq")
## [1] '1.32.0'
packageVersion("ggplot2")
## [1] '3.3.2'
data(enterotype)

set.seed(711L)

enterotype = subset_samples(enterotype, !is.na(Enterotype))

The Plot_Net Function

plot_net(enterotype, maxdist = 0.4, point_label = "Sample_ID")

Map as Color and Shape

plot_net(enterotype, maxdist = 0.3, color = "SeqTech", shape="Enterotype")

The Plot_Network Function

ig <- make_network(enterotype, max.dist=0.3)
plot_network(ig, enterotype)

Map Samples Variables as Color and Shape

plot_network(ig, enterotype, color="SeqTech", shape="Enterotype", line_weight=0.4, label=NULL)

Lower Maximum Distance

ig <- make_network(enterotype, max.dist=0.2)
plot_network(ig, enterotype, color="SeqTech", shape="Enterotype", line_weight=0.4, label=NULL)

Replace Jaccard with Bray-Curtis

ig <- make_network(enterotype, dist.fun="bray", max.dist=0.3)
plot_network(ig, enterotype, color="SeqTech", shape="Enterotype", line_weight=0.4, label=NULL)